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Abstract 

We present a theoretical agent-based model of cell evolution under the action of cytotoxic treatments, such as ra- 
dioteraphy or chemoteraphy. The major features of cell cycle and proliferation, cell damage and repair, and chemical 
diffusion are included. Cell evolution is based on a discrete Markov chain, with cells stepping along a sequence of 
discrete internal states from ’normal’ to ’inactive’. Probabilistic laws are introduced for each type of event a cell 
can undergo during its life cycle: duplication, arrest, apoptosis, senescence, damage, healing. We adjust the model 
parameters on a series of cell irradiation experiments, carried out in a clinical LINAC at 20 MV, in which the dam¬ 
age and repair kinetics of single- and double-strand breaks are followed. Two showcase applications of the model 
are then presented. In the first one, we reconstruct the cell survival curves from a number of published low- and 
high-dose irradiation experiments. We reobtain a very good description of the data without assuming the well-known 
linear-quadratic model, but instead including a variable DSB repair probability, which is found to spontaneously sat¬ 
urate with an exponential decay at increasingly high doses. As a second test, we attempt to simulate the two extreme 
possibilities of the so-called ’bystander’ effect in radiotherapy: the ’local’ effect versus a ’global’ effect, respectively 
activated by the short-range or long-range diffusion of some factor, presumably secreted by the irradiated cells. Even 
with an oversimplified simulation, we could demonstrate a sizeable difference in the proliferation rate of non-irradiated 
cells, the proliferation acceleration being much larger for the global than the local effect, for relatively small fractions 
of irradiated cells in the colony. 

Keywords: 


1. Introduction 

The development of cancer in a living organism fol¬ 
lows complex paths, not without seemingly contradic¬ 
tory features. From the broad point of view of sys¬ 
tems theory, the emergence of cancer cells appears 
as a stochastic process, in which one single cell sud¬ 
denly changes its nature without a traceable cause-effect 
mechanism, and starts proliferating abnormally. On the 
other hand, the rapid and uncontrolled development of 
tumoral tissues is unlikely to be described as a purely 
stochastic phenomenon, and metastatic propagation is 
drastically far from a simple diffusional flow: these fea¬ 
tures rather have the character of self-organising, non¬ 
equilibrium dynamical systems, susceptible of taking on 
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a chaotic or avalanche pattern under the effect of a small 
perturbation. 

Introducing the language of theoretical physics in 
the domain of cancer may seem unusual. However, 
physical-mathematical models have already accumu¬ 
lated a considerable tradition in cancer studies. Early 
analytical models based on coupled partial differential 
equations (see, e.g., |Brunton & Wheldon| ( p^80| ); [Sachs 


et al.| ( |20Q1| )), have been accompanied in recent years. 


and often superseded by complex numerical simulations 
models ( [Edelman et al.| |2010t [Deisboeck et al.| |2011t 
[Lowengrub et al. 2010[|Tracqui 2009| ), which attempt 
at following the space- and time-dependent dynamics of 
cancer growth, by adding an increasing wealth of de¬ 
tails and phenomenological correlations coming from 
biochemical and clinical studies. 

Despite the considerable efforts in modelling, cancer 
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treatments are still relying on a substantially empirical 
knowledge. Radiotherapy employs ionising radiation to 
eradicate cancer cells, mainly through the generation of 
DNA double-strand breaks (DSB), although the detailed 
mechanisms by which DSB and other sub-cellular le¬ 
sions are generated are still quite far from clear. Empir¬ 
ical descriptions, such as the so-called linear-quadratic 
model ( |Dale| 1985 ) are still extensively used in radio¬ 
therapy, to describe cell damage upon the delivery of 
ionising radiation, together with extensions, such as 
the ’Tumor Control Probability’ model ( [Kutcher 1996| ), 
aimed at predicting the clinical efficacy of radiothera- 
peutic protocols. However, a detailed correlation be¬ 
tween the radiation dose and its microscopic outcomes, 
both at the cell and tissue level, is still missing. 


In this work we develop, implement, calibrate, and 
apply a discrete-cell model with internal degrees of free¬ 
dom, describing both normal and abnormal cell evolu¬ 
tion, and accounting for localized damage and repair 
mechanisms, with the aim of studying the long-term 
evolution of a cell population subject to cytotoxic ther¬ 
apeutic treatments. Our main interest and application 
concerns the immediate and delayed action of radiother¬ 
apy treatments. However, the formalism developed here 
is enough general to be easily applicable to other cyto¬ 
toxic agents, such as chemotherapy, oxidative poison¬ 
ing, environmental (UV) radiation damage, and so on. 


The virtual cell population is represented by a large 
assembly (up to several millions) of individual stochas¬ 
tic age nts (see e.g. I Byme et ah] ( |2009| ); |Wang etaL 
( |2015| ); |Cilfone et aLp2015| )), endowed with a number 
of probabilistic properties (phenotypes), which allow to 
follow the evolution of individual cells through their 
daily cycle over very long time scales (days, months, up 
to years). Each cell has a local clock which goes through 
the Gl, S, G2 and M phases, typically (but not nec¬ 
essarily) following a 24h cycle. Cell duplication with 
inheritance is allowed, with individual probability laws 
depending on the cell state at any given time. Normal, 
stem, or tumor cells of various kinds can be included. In 
the simplest implementation, adapted to mimicking in 
vitro experiments on cell colonies, simulated cells live 
on a two-dimensional fixed square grid and can migrate 
by vicinal displacements. Diffusion of, e.g., oxygen, 
nutrients, or other chemical species is allowed on the 
same grid. 


These model cells can absorb a number of different 
damaging events, described by a Markov chain which 
changes the state of each cell from healthy, to pro¬ 
gressively damaged, to inactive. In this first paper we 
model only radiation-induced DNA damage in the form 


of single-strand and double-strand breaks. However, 
other DNA lesions, such as base excision, polymerisa¬ 
tion, clustered defects, could be included by extending 
the model, as well as damage to other vital cell com¬ 
ponents, such as mitochondria. Repair mechanisms are 
included aside of the damage, by assigning to each cell a 
set of probabilities to move its individual state upwards 
in the Markov chain state. Cells can also exit to a qui¬ 
escent, or a senescent state, which then follow special 
paths. 

Such a model will have predictive capability, after be¬ 
ing carefully calibrated on known evolution patterns of 
real biological cell lines. It should be able to predict 
the long term evolution (ranging from week-months, up 
to 5-10 years time, well-beyond the time scales accessi¬ 
ble to direct biological experimentation) of a population 
of cells with rather arbitrary characteristics, eventually 
subject to a series of time-defined treatments, for exam¬ 
ple by radiological or chemical agents, which may alter 
the cell vital cycle, notably by modifying its health con¬ 
dition and repair capabilities. In the present work, in 
order to calibrate the probability of inducing a SSB or 
a DSB, we performed photon-beam irradiation experi¬ 
ments in a clinical LINAC, on cultures of normal human 
dermal fibroblasts. The DNA damage was quantita¬ 
tively analyzed by searching foci of XRCCl and 53BP1, 
two proteins involved in the repair of SSBs and DSBs 
respectively, by means of immunofluorescence. 


While biophysical models covering some of the 
above characteristics have been already introduced 
in the literature, based on various mathematical ap- 
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proaches (see e.g. Sanchez-Reyes 


Stewart 


( |2001| ); [Kansal et al.| ( |2000| ); [Deisboeck et al.n2Qll| ); 
[Wang et ah ( |2015| ), and references therein), our model 
aims at assembling the most relevant features, in the 
attempt to develop a realistic platform for the virtual 
modelling of the long-term evolution of cell prolifera¬ 
tion and damage, following various types of therapeu¬ 
tic treatments. It is worth noting that the coupling of 
agent-based models with radio- or chemio-therapy, to 
simulate the action of external agents on cancer growth 
and/or arrest, is not yet fully developed. Even the most 
recent attempts in this direction (see e.g. |Kempf et~aL 


( |2013| ); [Powathil et al.| ( |2013| )) did not include an ex 
plicit simulation of the radiation damage, but rather as¬ 
sumed a pre-existing damage model (such as the linear- 
quadratic, etc.). An original contribution of the present 
simulation model is the introduction of explicit damage 
accumulation and repair, at the single-cell level. 

After an ample Section 2 describing the key details 
of the model, in Section 3 we include two benchmark 
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applications, aimed at demonstrating some key features 
of the computer model, namely: the reproduction of cell 
survival curves from irradiation experiments, and a sim¬ 
ulation of the so-called ’bystander’ effect. Note that, 
for the purpose of this first work, the two applications 
must be intended only as test cases, with no presumption 
of going in depth into the complex biophysical founda¬ 
tions, nor the medical implications of the corresponding 
phenomena. 


2. Model 


We use an agent-based model to describe the evo¬ 
lution of a cell population, under normal conditions, 
or in response to a physical or chemical perturbation. 
Each cell is represented as an ’agent’, endowed with 
probabilistic rules to follow in the course of the sim¬ 
ulation, both for their behaviour as independent units, 
and in interaction with each other. In our model, agents 
live on a fixed lattice (in the present study simply two- 
dimensional (2D) with fourfold symmetry), and can 
move on the lattice sites, carrying all the information 
about their state. A cell at a site i will be indicated 
as u{i). The time-dependent behaviour of each cell is 
characterised by a number of descriptors, collected in 
a state vector n(r) (see below). Such an implementa¬ 
tion is different from other lattice-based models, e.g. 
of Potts or lattice-gas type (for a review see [Anderson 


[]( |20Q7 1), in that the properties are dynamically carried 
by the agents, and not statically attributed to the lattice 
sites. We believe such a setting to be more fiexible in the 
numerical implementation, notably in view of the future 
developments of the model. 


Cells-agents may receive signals and input both from 
the environment and their neighbouring agents, as well 
as transmit signals to the environment and their neigh¬ 
bours. Agents make decisions based on both their in¬ 
puts and internal state. In future developments, cells 
will also include layers of subcellular decision-making 
rules; however, here we will only consider a constant 
set of rules, valid for the ensemble of cells. An agent 
may change state, proliferate, or undergo apoptosis or 
necrosis in response to surrounding conditions. Cellu¬ 
lar division requires space and a sufficient level of nutri¬ 
ents and oxygen. A cell may enter into a quiescent state, 
if the space is restricted (confluence), or if the nutrient 
supply is scarce. After cytotoxic treatments, a cells can 
arrest its life cycle. Quiescent, hypoxic or arrested cells 
can undergo apoptosis with a finite probability, after a 
defined length of time. Figure illustrates a fiowchart 
of the cell phenotype decision process. 


2.7. Normal cell cycle 

Aside of the global simulation time, hereafter indi¬ 
cated by t, each cell has a local clock = t — with 
^0 being the time of its last duplication, running over 
a cycle of 24h. The local time spans the four typi¬ 
cal phases of the cell cycle, namely: G1 from to 
750 min; S from td=15l to 1250; G2 from ^^=1251 
to 1370; M from ^^=1371 to 1440 min. (The timings 
are attributed conventionally.) Cells can be synchro¬ 
nised if needed, but the typical starting configuration is 
obtained with all cells randomly distributed among the 
four phases. Cell duplication should occur in the phase 
M, however a degree of randomness is allowed, by in¬ 
troducing a duplication probability: 

Pdupl = 4 (^j 4 .g-(«^-l 370 )/T 



with the duplication time constant T ^30 min, and (f) = 
1,2,3,4 for Gl, S, G2, M, respectively. Fluctuations in 
the duplication time of typical fibroblast cells, as well 
as for several other cell types including bacteria, are of¬ 
ten described by a ’’shifted-gamma” probability distri¬ 
bution (see e.g. Kutalik et~ar] (j 1 996|) ; [Rubinowj (j 1 996|) ; 


Stukalin et al. ( [1996 )), whose corresponding integrated 


cumulative distribution function is very close to a Fermi 
function. We use the latter in Eq.Q only because it 
is mathematically easier, however the numerical differ¬ 
ence with the shifted-gamma is practically irrelevant. 
The dependence on the parameter (j) is also rather ir¬ 
relevant numerically, since in all cases the probability 
function ^ is practically equal to zero, until is very 
close to 1370 min. (begin of M phase), around which 
time Pdupl rapidly goes to 1. 


The new (daughter) cell is spatially situated next to 
the old one, either in an empty lattice site, or by shift¬ 
ing nearby cells to empty sites in order to make room 
for the new one. When a cell duplicates, both its local 
clock and that of the daughter cell are reset to tQ=t. It 
is possible for a cell to arrest duplication by entering in 
the GO quiescent phase, for which 0=0. 


Normal cells may also enter a ’senescent’ state (dis¬ 
tinct from GO) some time after their development. 
Among other phenotype changes, senescence modifies 
or suppresses the cell duplication capability, while con¬ 
serving most of its metabolic activities ( |Hayf[ick[| 19651 
Cristofalo et al.|[2QQ^ , This can occur independently on 
the state of damage, in fact also normal, non-irradiated 
cells can naturally undergo senescence, under specific 
conditions [Refs Corinne] . This feature is taken into ac¬ 
count by counting the number of duplications for each 
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Figure 1: Flowchart of the cell decision algorithm and the Monte Carlo simulation. Black lettering indicates the actions of the computer code; 
red lettering indicates the corresponding biological/phenotype outcomes; in white lettering, the explanatory notes. The irradiation phase produces 
cell damage. In the subsequent evolution stage, while cells loop in the normal 24h cycle, they can repair damage, while following their doubling 
pattern. In this stage cell arrest can occur because of excess damage or senescence, followed by cell death; under given conditions restart can occur 
that, in turn, can yield back a normal cell, or evolve into a neoplastic cell. The coupled irradiation-evolution stages can be repeated cyclically, to 
simulate fractionated irradiation treatment. 


cell, Z), and introducing the maximum duplication num¬ 
ber Ds after which a cell enters into senescence, and the 
number Dq at which the population starts having senes¬ 
cent cells. Then, the probability P^upi is multiplied by a 
factor S=l for D <Dq, 5=0 for D > Ds, and: 


for Do < D < Ds. The values of the constants Dq and 
Ds are to be adjusted so as to reproduce typical senes¬ 
cence rates. Without loss of generality, in the following 
we will adopt Do=30 and Ds=60, a value assessed for 
the case of human skin fibroblasts [Refs Corinne]. The 
form of Eq.Q suggests that any cells that proliferated 
the maximum number of times Ds are arrested, therefore 
the ’senescent’ state may seem redundant with the GO. 


However, the GO can be entered at any time (for exam¬ 
ple because cells attained confiuence), and moreover the 
senescent state, notably in the case of keratinocytes, can 
be transient. Therefore, keeping this state allows cells 
to restart into a novel state even at (much) later times. 

2.2. Diffusion 

Chemical species /i with varying concentrations 
are allowed to diffuse on the lattice sites, and are accu¬ 
mulated at each time step Af in each cell u{i), instanta¬ 
neously located on the lattice site /, according to deter¬ 
ministic gradient fiow (Pick’s law): 
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where the sum runs only on the sites / nearest neigh¬ 
bours of site /, according to the chosen lattice topology. 
Since the neighbours realize in this first implementation 
a 2D diamond topology, the right-hand side of the pre¬ 
vious equation is just a discretized representation of the 
Laplacian, multiplied by the time step, with the diffu¬ 
sion time-scale 6^ playing the role (with appropriate 
dimensions) of a nominal diffusion coefficient. 

A source can be located at one specific lattice site, 
or an ensemble of sites, and diffuse through the empty 
lattice until reaching the cells. Or, it can be contained in 
a cell (for example, a secreted factor), diffuse on the lat¬ 
tice, and move along with the displacement of the cell. 

The cell membrane represents a semi-permeable bar¬ 
rier to almost all molecules and ions, with permeability 
coefficients much smaller than the diffusion in the sur¬ 
rounding fluid phase. Except special cases, therefore, 
diffusion on empty sites is considered instantaneous, 
while the membrane crossing of species jl is charac¬ 
terised by a diffusion time 6^ depending on the cell 
state and local density. 

It is worth noting that the local concentration of a 
species is accumulated in the cell u{i) occupying the site 
/, and not on the site itself. Concentrations can build up 
from lattice-diffuse sources as well as from cells. For 
example, oxygen concentration at cell u{i) is the result 
of the concentration field existing at site /, plus the even¬ 
tual gradient coming from the neighboring cells u{f). 
In the special case of oxygen diffusion, the concentra¬ 
tion at cells deep inside the colony can be considerably 
smaller than that of cells at the periphery, exposed to a 
much higher oxygen fiux. Similar examples could be 
the highly variable metabolic rate of cancer cells as a 
function of their distance from blood vases, or the com¬ 
plex chemistry of radiation-induced ’’bystander” effect; 
this latter will make for an application example in the 
following Sect. 3. 


2.3. Radiation damage 

Ionising radiation tracks deposit their energy in the 
cell in a time of ^10“^^ s. This energy density, in the 
form of secondary radiation, produces a variable ioni¬ 
sation density on its way, over a time of ~10“^ s, and 
over a length which, depending on the primary radia¬ 
tion nature and energy, can range from a few jimto sev¬ 
eral cm. The correlation between radiation energy de¬ 
livery and ionisation density is the linear energy trans¬ 
fer (LET). High-LET radiation, such as protons or alpha 
particles, have a high ionisation density along their path, 
however running along very straight tracks. Conversely, 


low-LET radiation, such as photons or low energy elec¬ 
trons, produce a much lower density of ionisation events 
per unit length, however their path is very diffuse and 
random, therefore their energy deposition is more ho¬ 
mogeneous in the cell volume. Anyway, only a rela¬ 
tively small fraction of the actual damage is produced 
by direct ionisation events, while the largest part is due 
to the secondary chemical species produced by the ra¬ 
diolysis of the molecules in the cytoplasm, i.e. mostly 
water, which thereby liberates free radicals OH* and H*, 
as well as free electrons (which go quickly into a sol¬ 
vated state). Such highly reactive species diffuse and 
attack the DNA {indirect damage), producing breaks in 
one (SSB) or both (DSB) the phosphate backbones, over 
a time scale of ^10“^-10“^ s. 


The ensemble of such events in our model is con¬ 
sidered instantaneous, and gives rise to stochastic dam¬ 
age events in the population of cells, according to rules 
which will be detailed in the following. The minimum 
time scale we consider is of seconds, for the irradia¬ 
tion time, and minutes for the follow up. Irradiation can 
take place according to different protocols, however it is 
usually delivered in batches of several tens of seconds at 
most. It is known that cell repair enzymes start working 
at considerably longer times (see e.g. 

( |2013| ); [Calini et al.] ( |2Q15| ), therefore 
that no cell repair activity takes place during irradiation. 


Georgescu et al. 
we will consider 


Radiation events are thought to follow a Poisson pro¬ 
cess ( |Kellerer||1996| ), therefore it may be justified to de¬ 
scribe the induced damage (both direct and indirect) as 
a Markov chain ( |Albright[ |1989t [Sachs et ak) |1989| ). 
We assume that cells can be in any state n G [0,m], 
with = 0 corresponding to a healthy cell with zero 
accumulated damage, and ^ = m to a cell with a max¬ 
imum of accumulated damages. The ^-th state of the 
cell is described by a set of indicators, or state vec¬ 
tor, n = for the v = 1, 

different types of lesions (in the present case, it is only 
V=1 for DSBs and v=2 for SSBs); i is the lattice site oc¬ 
cupied at time 0 is the cell phase; A is the cell state, an 
index 1,2,3,... for ’normal’, ’arrested’, ’dead’, ’neoplas¬ 
tic’ (plus additional labels for stem, cancer cells, and so 
on); Zv is the number of accumulated damages of type 
v; pv is the corresponding damage probability; ry the 
repair probability; c^ and 6^ the concentration and dif¬ 
fusion time of species jl. 

Let us define the probability Pn{t) that a cell is found 
in the state n at time t, with Pn>0 and = 1. The 
equation for a continuous-time Markov chain is: 


^^='ZWniPi{t) + S„ (4) 
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expressing the fact that the probability of observing the 
state n of the cell at time t is given by the sum over all 
possible probabilities of coming from any state / n, 
multiplied by the transition matrix Wni- Sn is a ’’source” 
term, representing the probability that a cell is arriving 
in state n from a cell division at time t. The transition 
matrix sums up all the processes leading to cell state 
evolution: 

= (5) 

V 

where represents the probability of making l-n=d 
lesions of type V, induced by a dose/unit time A (the 
dot indicating time derivative), and the probability 
of repairing d lesions of type V. 

When considering that the number of radiation events 
per unit time in each cell nucleus is relatively small, the 
matrix elements of D can be assigned the form of a Pois¬ 
son’s distribution. By assuming an energy spectrum of 
the radiation /(e), and by following the reasoning of 
[Sachs et ar] ( |2001| ), for each type of damage v the prob¬ 
ability of going from a state with n lesions to a state with 
n-\-d lesions is written in matrix form as: 



Zi, n. of DSB/cell 

Figure 2: Plot of the fraction of accumulated number of DSB le¬ 
sions/cell, zi, at different average dose levels from 0.5 to 2 Gy. Total 
number of cells, Nc=500. Symbols represent the raw results of the 
Monte Carlo simulation; continuous curves represent the fit with a 
Poisson probability law. 
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K,n-d = -rvd (9) 
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which will result in an exponential repair probability 
distribution. Of course, there would be little difficulty 

M. = j 

0 k\ 

f{e)de 


(7) 

in including more sophisticated models of repair kinet- 
ics (see, e.g., Cucinotta et al.lpOOS])). 


For numerical efficiency, however, the energy spec¬ 
trum is discretized, and for each energy interval the 
Poisson distribution is replaced by a limiting binomial 
( [Keinj et ari|2011| ): 


Lnn+d= 1 ™ 




( 8 ) 


with py = py{£) a piecewise, energy-dependent damage 
probability. By performing test simulations of irradia¬ 
tion at 2 Gy, we observe that already for mo ^ 5m the 
binomial is practically superposed to a Poisson curve 
for all values of dose (see Figure]^ showing the distri¬ 
butions of average DSB/cell as a function of the dose). 

The choice of the repair probability matrix R depends 
on the details of the enzymatic processes leading to 
strand rejoining after SSB, DSB, etc. In this first pa¬ 
per, we take a simplistic approach by choosing a simple 
linear model, in which the repair fraction is dependent 
on the number of lesions still present, with a generic re¬ 
pair probability r, whose numerical value is different for 


At this stage, our description of the relative biological 
efficiency (RBE) of the radiation is very simplified. To 
name a few important ingredients, we ignore complex 
damage types such as clustered defects, interactions be¬ 
tween defects, evolution of one defect type into another. 
However, such features can be added in further develop¬ 
ments of the numerical scheme; we are omitting them at 
present only to make for a more clear, albeit simplified, 
analysis of the results. 


The probabilities py and ry can be made dependent 
on the state of the cell. The damage probability can 
change according to the different radiosensitivity of a 
cell, e.g. increase because of an increased oxygen con¬ 
centration. The current version of the model includes 
the possibility of decreasing the damage probabilities 
Py by a simple sub-linear dependence on the total ac¬ 
cumulated dose. However, the radiation resistance is 
affected by many different competing factors, most no¬ 
tably a low oxygen level in the tumour mass ( |Graeber| 


et al. 2000). Therefore, a more complete description 


of such metabolic reactions should be accounted for by 
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means of appropriate parameters, controlled by a set of 
differential equations (see e.g. Wang et al. ( |2015 )). The 
repair probability may be made variable as well, for ex¬ 
ample the recruitment of repair proteins could be slowed 
down in a heavily damaged cell nucleus. Note that, for 
the sake of simplicity, in the examples of Sect. 3 we will 
only use constant damage and repair probabilities. 


In most experiments (see also Sect. 2.6 below) it is 
observed that similar cells, even belonging to the same 
line and culture, display variable radio sensitivity and re¬ 
pair capability. To mimic this effect, a specific depen¬ 
dence of the repair capability on the cell phase is in¬ 
cluded (see e.g. [Dikomey et al.| ( p^98| )), by setting: 


n 


(l)r if (j) <2 
Jr ii(l)>2 


( 10 ) 


for DSB repair kinetics, based on the fact that once the 
amount of DNA starts to increase in the cell nucleus be¬ 
fore mitosis, the cell repair machinery becomes quan¬ 
titatively less effective (for simplicity we take r 2 =const 
for SSBs, since they will not be relevant in the following 
simulation examples). It is worth noting that for can¬ 
cer cells of various types, the situation is more complex 
than for normal cells, in that a prevalence of radiosen¬ 
sitivity on the initial number of DSBs is more often ob¬ 
served, next to a reduced or even null dependence on the 
repair capability ( [El-Awady et al. |2003| ). 


2.4. Monte Carlo simulation 

The simulation space is defined by a 2D square lat¬ 
tice of N xN sites, uniquely labelled by a positive in¬ 
teger i <N^. In the present work we adopt a minimal 
(also called ’Von Neumann’) neighbourhood relation¬ 
ship, namely each site i interacts with the four neigh¬ 
bours (empty or occupied) located immediately above, 
below, left and right, in the 2D square topology. 

An initial number of cells Nc, typically quite smaller 
than the ensemble of lattice sites, is dispersed on the 
lattice, either at random, or in one or more compact 
colonies. Cells can move on the lattice by discrete 
jumps to neighbouring sites. It is already known from 
previous studies that the underlying lattice topology 
may induce unrealistic features in the tissue morphol¬ 
ogy. However, in the present work cell displacements 
are either not considered at all, or allowed only in order 
to attain confiuence; therefore this will not be a crucial 
limitation. On the other hand, if one is interested in sim¬ 
ulating tissue and tumour morphogenesis, a denser lat¬ 
tice topology or a Voronoi tessellation may be adopted 


( Moreira & Deutscfi||2002| ), and eventually a 3D exten¬ 
sion of the model becomes necessary. 


In the present version of the computer code imple¬ 
menting the agent-based model, different types of cells 
can be defined according to their phenotypes. Com¬ 
pared to normal cells, one could define stem cells, which 
follow a peculiar duplication pattern in that a stem cell 
divides into a normal and another stem cell; and tu¬ 
mour cells, which are defined according to a special 
set of the state vector n, e.g. an accelerated duplica¬ 
tion probability. In the examples to follow in Section 
3, we will only use either normal cells receiving a ho¬ 
mogeneous dose of radiation, or a mixture of irradiated 
vs. non-irradiated cells; two types of defects (DSB and 
SSB) will be tracked, although the cell death probability 
will be depending only on the number of accumulated 
DSBs. Such properties of the cell population can be 
easily modified, according to the subject of study. 


Solving analytically the Kolmogorov, forward-type 
equation 0 for a fully coupled set of space- and time- 
dependent probabilities is practically impossible. We 
therefore make recourse to a Monte Carlo stochastic 
method of solution (see Figure [^. At each discrete up¬ 
date by Ar of the global time t, all the cells are scanned, 
and their state vector n is updated. Probabilities for 
the various events a cell may undergo are sampled ac¬ 
cording to a rejection technique, i.e., a random number 
^ G (0,1) is drawn from a fiat probability distribution, 
and compared to the event probability P. The event is 
accepted if <^ < P, otherwise it is rejected. For numer¬ 
ical efficiency, all such events are described in the fol¬ 
lowing by continuous probability functions, however it 
may be noted that the functional shape practically corre¬ 
sponds to binary, ”yes/no” options, at each Monte Carlo 
sampling step. 


For example, let us focus on a particular cell i e Nc'. 
as long as the time t increases, its local clock t^ ad¬ 
vances, and its duplication probability P^upi increases 
from nearly 0 in the G1 phase, to nearly 1 in the M 
phase. Correspondingly, at each time step t, the new 
0< <1 will stochastically sample the increasing prob¬ 

ability Pduph by producing a duplication event at a ran¬ 
dom time ^0 distributed according to Pduph The same 
random sampling happens for the probability of going 
from Zv to Zv lesions (with v=l for DSBs and v=2 
for SSBs, and damage probability py ifd> 0, or repair 
probability ry if J < 0); for the probability of going into, 
or escaping from, the quiescent state GO; and so on. 


In a typical simulation of an irradiation cycle, cells 
are irradiated for some time of the order of 10s to 100s 
of seconds, with a time-step A^=l s. Then, the evolution 
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of the partly damaged population is followed for a time 
of the order of several hours, with a typical Ar=60 s or 
larger. Eventually, the irradiation-evolution cycle may 
be repeated, to simulate fractionated radiotherapy. For 
the sake of simplicity, we assume that no cell evolution 
takes place during the short irradiation time, however 
this limitation can be easily removed. 

As a result of radiation induced damage, cells can 
arrest the normal duplication cycle in any of the four 
phases, because some cycle checkpoint is not com¬ 
pleted. We construct for this event a probability: 


1 if zi > 0 
l-S ifzi=0 


( 11 ) 


meaning that any cell having undergone DSB damage, 
or any cell entering senescence (S > 0, see Eq. 0) can 
be arrested. 

Once in this arrested condition, a cell can take at each 
time step one of three alternatives (Fig. [^. Firstly, it 
can die and be eliminated by apoptosis: for example, 
it is known that the p53 protein, usually a regulatory 
transcription factor, can also trigger the expression of 
genes inducing apoptosis, when excessive DNA damage 
accumulates ( |Roos & Kaina|[2006| ). For this event, we 
introduce a probability: 


Pdeath = min 



( 12 ) 


with a > 1 a ’retarding’ parameter, and Ncrit a criti¬ 
cal threshold of accumulated DSBs leading to apoptosis 
(’mortality’ parameter). The value of the threshold de¬ 
pends on the cell type and on the irradiation conditions, 
typical values being in the range Ncrit ^ 5 ^ 50. 

Alternatively the cell can restart its life cycle, with 
probability: 

Prestart = (13) 

increasing at a rate l/j8 ^ 0.5 to 2, expressing the effi¬ 
ciency of the repair action, and getting equal to 1 when 
Z\ is back to 0. Or, finally, the cell can remain in the qui¬ 
escent state with probability (1 — Pdeath — Prestart) until 
the next time step. 

However, it is known that the reentry in the life cy¬ 
cle of a damaged cell can lead to a potentially cancer¬ 
ous cell, especially at lower levels of damage and at¬ 
tack to the tumour suppressor genes, events not imme¬ 
diately leading to a tumoral cell, but rather inducing a 
genomic instability, which is only one of the required 
components likely to induce a cancer. Such cells may 
be indicated as ’neoplastic’. 


2.5. Calibration by irradiation experiments on human 
fibroblasts 

Normal human dermal fibroblasts (NHDFs) used in 
this study were PromoCell F-IMC from a 1-year-old 
Caucasian male. Cells were cultured in incubator at 
5% CO 2 and 37 °C, in basal medium (FBM - Fibrob¬ 
last cell Basal Medium, Fonza) with 2% FBS (Foetal 
Bovine Serum), plus human fibroblast growth factor 
(hFGF), insulin at 5 mg/ml, and antibiotics Gentam¬ 
icin 50 /ig/ml and Amphotericin B 50 /ig/ml. Cells 
were plated at 200,000 cells per 100-mm Petri dish, and 
cultures were always split at 70% confluence with the 
Reagent PackTM by Clonetics (Hepes buffered saline 
solution, Trypsine/EDTA 25 mg and 10 mg per 100ml, 
TNS Trypsin Neutralizing Solution). Population dou¬ 
bling was calculated at each passage, after cell counting 
with a Thoma or Malassez counting chamber. The num¬ 
ber of population doublings, PD, was calculated as: 

pj^ ^ collplat) 

ln2 ^ ^ 

^colh ^plat being the number of collected and plated 
cells, respectively. 

Irradiations took place at the Varian Primus CFINAC 
of the ’’Oscar Fambret” center. Cells were cultured 
either in 96-well plate, deposited on 4 cm of a plas¬ 
tic plate (equivalent tissue) to ensure electronic equilib¬ 
rium. (Note that for a 20 MV accelerating tension, the 
photon spectrum is peaked at a much lower energy of 
about 3 MeV, with an average energy of about 10 MeV.) 
The photon beam was directed from below the table, 
with intensity adjusted to provide 200 monitor units at 
the isocenter of the culture dish. 

After the irradiation, cells were fixed with formalin, 
and permeabilised with triton 1% in PBS (Phosphate 
salt buffer: 3 . 2 mMNa 2 HP 04 , 0.5mM KH 2 PO 4 ,1.3mM 
KCl, 135mM NaCl, at pH=7.4). After washing, non¬ 
specific sites were blocked for Ih in 5% skim milk di¬ 
luted in PBS. Cells were then incubated in a solution 
of primary antibodies diluted at 1/100 or 1/200 (respec¬ 
tively for anti-53BPl and anti-XRCCl) in the buffer 
(anti-53BPl from Santa Cruz Biotechnology, SC-22760 
and anti-XRCCl from Santa Cruz Biotechnology, SC- 
11429). After several washes in PBS, cover slips were 
incubated in a solution of secondary antibodies A21206 
anti-IG (rabbit) from Fife Technologies, diluted at 1/500 
and coupled to a fiuorochrome (AFEXA FFUOR 488). 
Subsequently, cell nuclei were coloured by a solution 
of Hoechst 33342 (Sigma) 5 mg/1 in PBS. 96-well plate 
were finally analysed by HCS Operetta (Perkin Elmer). 
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Counting of XRCCl and 53BP1 foci was carried out 
with an image analysis software (Columbus by Zeiss). 


3. Test cases in modelling cell culture irradiation 


In this section we present two applications of our 
model, exploring only some of the wide range of pos¬ 
sible simulations of cell evolution that it enables. We 
stress once more that such test cases are not meant to 
explore the complexity of the corresponding biophysi¬ 
cal problems, but they only serve as the first, necessary 
test phase of the model. 


To begin with, we need to firstly adjust the model 
kinetics on our experimental results of cell irradiation. 
Figure shows with square symbols the results of the 
counting of DSBs (blue) and SSBs (red) after a nom¬ 
inal 2 Gy irradiation. Note that the experimental con¬ 
trol value (indicated by crosses at r=0) is practically 0 
for DSBs, while it is around 10 SSB/cell; such a back¬ 
ground corresponds to daily generation of damage from 
sources other than the irradiation. After subtracting the 
background, it is seen that over the 24 hours about 25% 
of the DSBs are still unrepaired (blue squares), while 
practically all the SSB have been repaired (red squares; 
white squares, before subtraction of control). However, 
it is worth noting that for SSB there remains a disparity 
between cells, as shown by the larger standard deviation 
in the histogram, some cells still showing ^25 XRCCl 
foci after 24h. The average value of about 6 DSB/Gy 
observed here is quite smaller than the accepted val¬ 
ues, in the range of 15-35 DSB/Gy ([Lobrich et al.[|1993 


Schultz et al.||2000l|El-Awady et al.[|2003| l. The source 

of such a discrepancy is not yet clear, it might be related 
to the particular kind of signalisation protein chosen for 
damage identification. 


For the computer simulation we will use a small test 
sample, made of Nc=500 cells distributed on a square 
lattice, and subject to a homogeneously distributed ra¬ 
diation dose. From our preliminary tests, such a size of 
population is already sufficient to obtain a good statis¬ 
tics on the system response. In the same Fig. we 
show two calculated histograms for creation (t < 0) and 
reparation (t > 0) of DSB (blue) and SSB (red). We cal¬ 
ibrated the probabilities pi and p 2 of generating a DSB 
or a SSB, respectively, as well as the corresponding re¬ 
pair probabilities ri and r 2 , by using the experimental 
data as a reference. For the sake of simplicity, we as¬ 
sume the dose of 2 Gy, at a rate of 3 Gy/min (lower 
range of a typical LINAC machine), to be delivered by 
a monochromatic spectrum of gamma-rays at 20 MV. 
(Note that the x-axis in the figure give the time in sec¬ 



Figure 3: Kinetics of the experimentally observed 53BP1 and XRCCl 
foci after 2 Gy irradiation (symbols), with repair probabilities from 
Eq. (histograms) adjusted to match the experimental kinet¬ 
ics. Crosses: background signal in absence of radiation (red=SSB, 
blue=DSB). White squares: XRCCl raw data; red: XRCCl-n after 
subtracting the background at t=0; blue: 53BP1 data. Note the larger 
error bars for SSBs, pointing to a larger dispersion of the cell popula¬ 
tion for this type of lesions. The blue and red histograms correspond 
to the simulated model kinetics, with damage probabilities {pi^pi}, 
and repair probabilities {ri, r 2 } adjusted to reproduce experimental ir¬ 
radiation data. The x-axis gives the time t in seconds, for r < 0, and 
in minutes, for f > 0. 


onds, for r < 0, and in minutes for ^ > 0). The best 
reproduction of the experimental data is obtained with 
Pi = 3.0 X 10-4, /72 = 6.5 X 10-4, n = 1.5 X 10-^r2 = 
2.5 X 10-^. The resulting histograms have a nearly ex¬ 
ponential decay (as expected from Eq.([^), with the re¬ 
spective f=0 values giving the pv, and the decay con¬ 
stants giving the ry. Also, note that since we have 3-f3 
experimental data points, there is no over fitting of the 
four parameters. 

Figure]^ shows the effect of using a repair probability 
that depends on the cell phase. The open squares for the 
’mix’ case, in which the repair probability follows Eq. 
( p^ , correspond to the best fit already given in the pre¬ 
vious Fig. For the other cases, it can be seen that by 
restricting the repair to S-phase only a slower recovery 
of damage is obtained, while, at the other extreme, the 
faster recovery would correspond to concentrating the 
repair only in the G2 and M phases. In the remainder of 
this work, we will use the ’mix’ option corresponding 
to the best fit. 

After calibrating in this way the kinetics of our 
model, we will present a first example in which we ap¬ 
ply the model to reproduce the well-known results of 
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Figure 4: Effect of including a cell-phase-dependent repair probabil¬ 
ity (see Eq.flO)). ’Mix’ is the average, with repair functioning in all 
cell phases, ’S’, ’Gl’, ’G2-M’ refer to a repair probability restricted 
only to the corresponding cell phase. Continuous lines represent ex¬ 
ponential fits to each set of symbols. 


typical cell survival curves, following irradiation with 
increasing levels of dose. Then, in a second example 
we will test the ability of the model to incorporate spa¬ 
tial information, by developing a set of simulations of 
the so-called ’’bystander” effect, namely the alteration in 
the response of non-irradiated cells when some nearby 
cells are irradiated. 


3.1. Survival curves for irradiated cells 

The ’’linear-quadratic” (LQ) model is the most widely 
accepted mechanistic model of cell killing by radia¬ 
tion (see e.g. |Dale| ( |1985| )). While being a strictly 
phenomenological representation of experimental data 
about cell survival as a function of the total dose, the 
LQ is often justified by the concept of binary DSB mis- 
repair: since it is believed that the main chromosome 
aberrations (’’dicentric”) should likely occur after the 
wrong rejoining of a pair of close-by DSBs ( [Wlodek 
& Hittelnian| |1988| ), it may occur that prolonged expo¬ 
sure to radiation would allow only one of the DSB to 
be repaired, before the second is generated. Although 
this is the most usual way to provide a motivation to 
the standard LQ approach, alternative biological expla¬ 
nations have also been advanced, such as the repair- 
misrepair ( [Tobias | [1985] ), the lethal-potentially lethal 
( jCurtisj p"9^ l, or the two-lesion kinetic model ( |Stew-| 
ait| 2001j ), which give practically the same numerical 
results of the LQ; as well as alternative models, vari¬ 
ously based on concepts of saturation of the repair ca¬ 
pacity ( jSanchez-Reyesj |1992t jCucinotta et al. |2008| ). 


Figure 5: Plot of the cell death fraction from a colony of iV=500 cells, 
as a function of the DSB repair probability ri, for simulated irradi¬ 
ations at dose levels between 1 and 10 Gy (indicated for each set of 
data), at a rate of 3 Gy/min. Damage probabilities pi,P 2 , and repair 
probability r 2 , are fixed at the same values reproducing experimental 
irradiation data in Fig. [^ Cells are observed for 24h after the irra¬ 
diation, and the mortality threshold is fixed at Ncrit=^ (see Eq. Gl)- 
Continuous curves are best fits with a sigmoidal function, for each set 
of simulation data. 


Overall, it is generally acknowledged that the biophys¬ 
ical basis of all such models rests on quite speculative 
assumptions about the microscopic events leading to the 
observed cell response. Their most important value, no¬ 
tably for the LQ model, rather resides in their empiri¬ 
cal utility in dose treatment planning (for a comparison 
among the various models, see e.g. [Brennerj ( |2008| )). 
However, one cannot say to have learned much about 
the biophysics of DNA damage from empirical models 
of this kind. 

We used our simulation model to attempt an inter¬ 
pretation of typical cell survival curves, however with¬ 
out including any of the above detailed (but specula¬ 
tive) microscopic assumptions about the typology and 
evolution of DSB generation. As an example, we will 
test the hypothesis ( jSanchez-Reyesj |1992t jCucinotta"^ 
ak] [2008] ) that the repair probability may be dependent 
on the dose, by means of some ’saturation’ mechanism, 
which could be further traced back to a molecular ori¬ 
gin. However, in our simulations we will not explicitly 
introduce the saturation hypothesis. Rather, we will fix 
some damage threshold Merit, and explore the cell sur¬ 
vival fraction as a function of the radiation dose and 
the DSB repair probability ri. The other ingredients are 
those already given in Sect. 2, namely DSBs are taken 
to occur independently, with Poisson-like statistics, fol¬ 
lowed by a simply exponential repair capability. No fur- 
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ther fitting is needed, since we use the values of prob¬ 
ability corresponding to the previous fitting 

(Sect. 3.1 above) of our own experimental data at the 
dose of 2 Gy. In this way, we aim at showing (i) that 
the model is data-independent (since we use data from a 
particular fitting, to reproduce data obtained in entirely 
different experiments), and most importantly (ii) that 
the model can autonomously generate a non-predefined 
response, thus providing a biophysical interpretation to 
the experimental data. 

Figure shows the evolution of the cell death frac¬ 
tion as a function of the value of ri, for irradiations 
at different dose levels A=l,2,...10 Gy, at the rate of 3 
Gy/min. The mortality threshold is fixed at a quite low 
value, Ncrit=^, meaning that on average 5 DSBs in a 
cell are already enough to lead to apoptosis. Cells are 
observed for a time of 24h after the irradiation. For the 
sake of simplicity, we are using a phase-independent re¬ 
pair rate, i.e. the ’mix’ in Fig. All the curves for the 
various simulations display the same behavior, namely 
the death probability is close to 1 for low values of n , 
and it rapidly decreases around a critical value of ri, de¬ 
pending on the dose. Upon increasing ri all the curves 
end up to zero mortality, meaning that the repair proba¬ 
bility gets large enough to provide a full healing of DSB 
damage within the 24h time window. Note that for the 
lowest dose simulated of 1 Gy, the plot does not start 
at 100, since the number of DSBs generated by such a 
small dose is below Ncrit=^- 


Now, we take some typical cell survival curves from 
published irradiation experiments, and compare the ex¬ 
perimental data to our simulated survival fractions from 
Fig. Namely, for each value of dose A, we extract 
from the corresponding curve in Fig. |^the value of ri 
which allows to best fit the experimental survival frac¬ 
tion at that dose. Figureshows three such experimen¬ 
tal curves. Two are from the work of |S0rensen et al.| 
( 2011| ), obtained by irradiation of V79 (blue) and FaDu 
(red) Chinese hamster cells with a LINAC, at dose rates 
of 5.01, 9.99 and 29.91 Gy/min (see left y-axis). The 
three colour curves represent our model fit to the exper¬ 
imental data, as proposed in the original works. In such 
experiments, no dependence of the survival fraction on 
the dose rate was evidenced, in fact the LQ fits for the 
three different rates are practically identical. Note that 
such dose rates are much larger than those used in clas¬ 
sical survival curve experiments by smaller radioactive 
sources: as an example we include in the Figure also 
data by Ma et al. ( |2013| ) (green) on XI Chinese hamster 
cells irradiated by a 250 kV x-ray tube. 

In the same Figureright axis, we also plot the val¬ 



dose (Gy) 


Figure 6: Left y-axis: plot of the experimental data of survival frac¬ 
tion for V79 (blue dots) and FaDu (red dots) cells, fro m|S0rensen et| 
[ar]j2011h and for XI cells (green dots) from |Ma et al.H2013) . Con- 
tinuous curves represent a fit to our simulation data. Right y-axis: 
square symbols: values of ri giving the best fits (color curves) to the 
experimental survival fraction, for each dose value; symbol colors cor¬ 
respond to each experimental and fitted curve. The dashed line is an 
exponential fit to the ri values from the model. 


ues of ri which provide the closest fit to each curve (blue 
squares for the V79, red squares for the FaDu, green 
squares for the XI), at each value of dose. It can be 
seen that the best values of ri are not at all constant, 
but follow a clearly saturating pattern, i.e. increasing 
steadily at low doses and aiming at a constant value 
(about ri ^ 2.4 x 10“^) at higher doses. The dashed 
curve in the figure is an exponential fit to these ri val¬ 
ues, of the type a — bexp{—cri), with <2=2.4, b=3 and 
c=0.25. 

The important point here is that our dynamical sim¬ 
ulation model is capable of producing survival curves 
in very good agreement with various sets of experimen¬ 
tal data, without making a priori assumptions about the 
biophysical mechanisms underlying. In a sense, we let 
the cells ’choose’ their behavior, and we look which par¬ 
ticular behavior corresponds to the observed experimen¬ 
tal data. In this way, the model suggests a biophysical 
mechanism underlying the observed shape of the sur¬ 
vival curves: rather than by making hypotheses of, e.g., 
linear-plus-quadratic components (of doubtful physical 
origin), the model only includes biologically observable 
damage (SSBs and DSBs, although only the latter are 
considered as lethal, according to the current knowl¬ 
edge). With this initial hypothesis, the model is able to 
suggest that cells may adopt a molecular repair response 
that is exponentially saturable. One could then further 
speculate about the molecular connection, for example 
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imagining that there is a finite supply of repair proteins 
in the cell nucleus, which fails when the damage is too 
extended. While we used here this set of computer sim¬ 
ulations only as an example of application of the model, 
a more complete study of this effect is underway, and 
we defer further considerations to a future work. 


3.2. Computational model for 'bystander' effect in ra¬ 
diobiology : long-range short-range diffusion. 


The so-called ’’bystander” action is seen to occur in 
experiments performed under widely different condi¬ 
tions, both in vitro ad in vivo, and it appears to fall under 
two broad categories: (i) local action over short cell-cell 
distances; (ii) global action over longer distances. 


The cellular effects of such bystander actions gener¬ 
ally involve an upregulation of the metabolism of co¬ 
affected cells, such as the oxidative pathways ( |Narayan| 
et al.[ \\991\ |Mothersill & Seymour 1998| ), the p53 
damage-response pathway ( |Azzam et al. |1998|), in - 
creased levels of interleukin-8 ( Narayan et al. 1999| ), 
AP endonucleas e (|Iyer & Lehn^ |2002 ), TGF j8l 
( Iyer & Lehnert| 2000[ ), and others. Notably, increase 
in cell proliferation activity has often been associated 


with such increased metabolism (Iyer & Lehnert 2000 


|2002[ [Gerashchenko & Howell[ [2003 [2005 j ). This latter 

finding motivates a second application of our model, in 
order to show the importance of coupling the damage- 
repair capability with diffusion and cell topology. 


We performed a series of coupled irradiation- 
diffusion simulations on a sample colony of A/c=500 
cells, with the same probability parameters as above. 
To simulate co-culture of irradiated and non-irradiated 
cells, randomly picked sub-populations corresponding 
to 5, 10, 20, 50, 75 and 90% of the colony were irra¬ 
diated with the same protocol of Sect. 3.1. The by¬ 
stander effect was simulated in the two variants above 
described, namely (i) local and (ii) global, by allowing 
the irradiated cells to diffuse a generic pro-mitogenic 
factor (such as TGF j3l), here labelled B for ’by¬ 
stander’ , whose concentration (P is supposed to increase 
over time in each cell, according to the diffusion equa¬ 
tion The source term was set to = 1 in irradiated 
cells, and = 0 elsewhere. The increase in prolifer¬ 
ation was linked to the concentration, by shifting each 
local cell clock as: 


id'— id 3-At{lyc^) (15) 

with 7 an empirical efficiency factor, ranging from 0.1 
to 2. In this way, cells that accumulate larger concen¬ 
trations of the ’B’ factor may anticipate their time to 


the next division, up to twice the normal rate. In the 
same spirit, we also tested the effect of a cytotoxic fac¬ 
tor (e.g., reactive oxygen species, hydrogen peroxide, 
nitric oxyde) leading to proliferation arrest, by setting 
7=-0.5, -1. Moreover, we tested the possibility that the 
’B’ factor may also induce apoptosis when exceeding 
some threshold (for example, c^ > 0.9), by relating the 
concentration to Pdeath Eq.(p^ as: 


Pd eat h 


1 

l^^-100(c^-0.9) 


(16) 


Then, the ’local’ variant of the bystander effect 
(i), was represented by assigning diffusion coefficients 
D^=\ and (in arbitrary units) to the irradiated 

and non-irradiated cells, respectively. The ’global’ by¬ 
stander effect (ii), was instead represented by assigning 
D^=\ to irradiated, and increasing values of 
10“^, or 10“^ to the surrounding, non-irradiated cells. 
An example of the typical effect of such choices is rep¬ 
resented in the 2D cell maps of Figure [7] 

The results of the simulations are summarized in Fig¬ 
ures Hand 121 The first one shows the duplication rate 
of non-irradiated cells after 24h, in a colony in which 
variable proportions of irradiated cells (in abscissa) are 
randomly mixed. The four panels correspond to the four 
values of :^0, =0.001, 0.01, 0.1. The curves in each 
panel correspond to different values of y, the coefficient 
relating the concentration of the B factor to the prolifer¬ 
ation rate. From such data, only minor differences be¬ 
tween the ’local’ and ’global’ bystander effect could be 
inferred: all the families of curves resemble each other, 
all showing a sizeable increase in the duplication rate 
as a function of y (the normal proliferation rate in the 
absence of irradiation would be 2 after 24h). For y ap¬ 
proaching zero, the proliferation rate goes back to con¬ 
stant, around its normal value of 2. Negative values of 
7 give instead reduced proliferation rates, with the cell 
population being reduced to some constant value, since 
part of the cells go into the arrested condition. The case 
in which the B factor induces apoptosis is the curve la¬ 
belled ’m’ (for ’mortality’) in each panel. It can be seen 
that in this case the cell colony can be completely killed 
after 24h. In this latter case, some correlation between 
the fraction of irradiated cells and the diffusion coeffi¬ 
cient exists: at the smaller values of , the mortality 
rate is nearly exponential with the fraction of irradiated 
cells, while at larger diffusion coefficients the mortality 
grows very quickly, and already a fraction of 10 to 20% 
of irradiated cells is sufficient to propagate very effec¬ 
tively the simulated cytotoxic factor. However, by com¬ 
paring the upper-left and the lower-right panel, at least 
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Figure 7: Concentration maps from the simulation of bystander effect in a colony of 500 cells with a randomly distributed 5% fraction of irradiated 
cells. Rows from top to bottom: diffusion coefficient for non-irradiated cells Z)^=~0, 0.001, 0.01, 0.1. The top row represents a local bystander 
effect, requiring direct cell-cell contact; the lower rows correspond to the global effect, with increasing extracellular diffusion velocities. Columns 
from left to right: simulation time r=0, 1, 5, 24h. The continuous color coding represents the concentration of the secreted factor in each cell (blue 
c^=0, white c^=0.5,red c^=l). 


a qualitative difference can be appreciated, the local ef¬ 
fect showing a dependence of the proliferation rate on 
the fraction of irradiated cells, which is practically ab¬ 
sent in the global effect at the highest values of . 

One has to look at Figure to obtain a much clearer 
evidence of a sizeable difference between the local and 
global versions of the bystander effect. The plot repre¬ 
sents two families of curves, each one corresponding to 
different values of y as in the previous Fig. How¬ 
ever, in this case the abscissa represents the values of 
in logarithmic scale. The two families of curves are 
for two extreme values of fractions of irradiated cells in 
the colony, a small 5% (red lines, dashed) and a rather 
large 75% (black lines, full). It can be clearly appreci¬ 
ated that, for a high fraction of irradiated cells, both the 
local and global effect give results that are practically 
independent on the value of the main parameter to 


determine the proliferation rate being only the value of 
y. However, at the smallest fractions, and for the larger 
values of 7 > 1, there is a marked dependence of the pro¬ 
liferation on T)^, the ’global’ bystander effect leading to 
much higher proliferation rate than the ’local’, indeed 
more than doubled for 7=2. This dependence seems less 
apparent for the negative values of y. 

Such a result of the model is not at all surprising, and 
could have been readily predicted simply on the basis 
of logical reasoning: if the number of ’source’ (i.e., ir¬ 
radiated) cells is large, there is practically no difference 
between a local and a global action; if the number of 
source cells is quite sparse, instead, a global action is 
necessarily more effective than a local one, for a given 
amount of time. Being only a test problem, this is an in¬ 
direct proof that the model works correctly. Clearly, our 
naive, two-parameter representation of the bystander ef- 
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Figure 8 : Plot of the duplication rate of non-irradiated cells after 24h, 
in a colony of 500 cells including fractions of irradiated cells ranging 
from 5 to 90%. In each simulation, the diffusion coefficient of irradi¬ 
ated cells is always Z)^=l. The diffusion coefficient of non-irradiated 
cells varies in each panel. Top row, left ~0, right =0.001; bot¬ 
tom row, left =0.01, right =0.1. In each panel, the various curves 
correspond to different values of 7 in Eq. 0 . 


feet is far from capturing the complexity of this elusive 
phenomenon, which requires instead a much deeper in¬ 
vestigation. We used it here, just for the sake of demon¬ 
strating the power of a theoretical model of cell evolu¬ 
tion capable of coupling damage and repair information 
at the single-cell level, together with spatial and topo¬ 
logical information at the multicellular level. 


4. Discussion and Conclusions 

In this work we developed and tested an agent-based 
model of cell evolution under the action of cytotoxic 
treatments, aimed at including the major features of cell 
cycle and proliferation, cell damage and repair, chemi¬ 
cal diffusion. In this first implementation of the model, 
the cells-agents live in a two-dimensional lattice, with a 
diamond-neighbour topology. Cell evolution is mathe¬ 
matically driven by a Markov chain, with cells taking on 
a series of discrete internal states from ’normal’ to ’in¬ 
active’ (or ’dead’). Probabilistic laws are introduced for 
each type of event a cell can undergo during its life cy¬ 
cle: duplication, arrest, apoptosis, senescence, damage, 
healing. The system is simulated with a time-forward 
explicit algorithm, with a Monte Carlo sampling of the 
various probability distributions. 

We firstly calibrated some of the free parameters with 



Figure 9: Plot of the duplication rate after 24h, in a colony of 500 cells 
including 5% (red, dashed) or 75% (black, full curves) fraction of irra¬ 
diated cells, as a function of the diffusion coefficient of non-irradiated 
cells. The diffusion coefficient of irradiated cells is always D^=l. The 
various curves correspond to different values of yin Eq. tn}. 


a series of cell irradiation experiments, carried out in 
a clinical LINAC at 20 MV and quite high fiuence 
(^3 Gy/min). Single- (SSB) and double-strand breaks 
(DSB) in irradiated cells were counted by means of au¬ 
tomated image analysis, after staining with XRCCl (for 
SSBs) and 53BP1 (for DSBs) proteins. The resulting 
damage and repair kinetics were deduced, and included 
in the model parametrization. 


The very first test of the model was directed at re¬ 
obtaining the cell survival curves from a number of 
published low- and high-dose irradiation experiments. 
Such curves are generally interpreted on the basis of the 
so-called ’’linear-quadratic” model, which presumes the 
DSB damage to occur in two different forms, a single¬ 
hit and a double-hit mode (hence the two terms in the 
LQ model). We could obtain a very good representation 
of the same cell survival curves, however without as¬ 
suming any special hypothesis about the types of DSBs 
or other. We only let the cell DSB-repair probability as 
a variable, and we obtained that in order to fit the ex¬ 
periments such a probability must saturate as a function 
of the increasing dose, as previously suggested e.g. by 


Sanchez-Reyes (1992) and Cucinotta et al. (2008). This 


simple test already shows the power of the simulation 
method, in that it is capable of suggesting a biophys¬ 
ical behavior that cells may adopt (an exponentially- 
saturating DSBs repair capability), instead of relying on 
various pre-conceived hypotheses. 


As a second test, we attempted to simulate the so- 
called ’bystander’ effect in radiotherapy, namely the 
possibility that non-irradiated cells located in proximity 
of irradiated cells may develop a sequence of metabolic 
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response similar to these latter. We tested the two 
extreme, opposing hypotheses of a ’local’ bystander 
effect, activated by short-range intercellular diffusion, 
versus a ’global’ effect, activated by the long-range dif¬ 
fusion, e.g. of some factor secreted from the irradi¬ 
ated cells in the extracellular fluid. Even if our sim¬ 
ulation was exceedingly simple, and many important 
factors were explicitly neglected, we could demonstrate 
some sizeable difference in the proliferation rate of non- 
irradiated cells, the increase in proliferation being much 
larger for the global than for the local effect, at least 
for a relatively small initial fraction of irradiated cells 
in contact with normal cells. 

In this first study, primarily devoted to introducing 
and testing our agent-based model, we did not use 
the information about the other type of lesions already 
included in the algorithm, namely the single-strand 
breaks, for which we anyway performed the experimen¬ 
tal calibration in the same way as for the double-strand 
breaks. We are however planning to use such a feature 
in a forthcoming study of tumour dormancy, as well as 
the other features not yet fully developed, such as the 
cell motility, which would be essential for example in a 
context of tumour heterogeneity studies. 
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Symbol 

Variable 

Units 

Role 

Values 

N 

lattice size 


linear size of the square lattice 

— 

i 

site index 


label of lattice site 


Nc 

number of cells 


running size of the cell population 

— 

u{i) 

cell index 


label of cell w on a lattice site 

[IM 

A 

dose 


amount of energy supplied over a given time interval 

[0,oo] 

n 

cell state vector 


vector containing all local cell parameters 

— 

t 

global time 

seconds 

universal ’laboratory’ time of the simulation 

1^ — oo^ coj 

At 

time step 

seconds, 

minutes 

discrete simulation time increment 

[0.001-0.1] 

td 

cell time 

minutes 

time clock local to each cell 

[0,1440] 

to 

duplication time 

minutes 

time since last duplication for each cell 

[0,00] 

T 

duplication time constant 

minutes 

see Eq.|T] 

30 

0 

cell phase index 


defines me cell phase (G0,G1,S,G2,M) 

[0,4] 

A 

cell state index 


normal, stem, cancer, arrested, dead, neoplastic 

[1,6] 

V 

damage type 


double- or single-strand breaks (DSB, SSB) 

[1,2] 

Zv 

damage counter 


counts accumulated number of defects of type v 

[0,oo] 

Pv 

damage probability 


probability of generating a defect of type v 

[0,1] 

ry 

repair probability 


probability of healing a defect of type v 

[0,1] 


critical damage 


number of lethal damage (DSBs) above which a cell undergoes 
apoptosis 

[5,100] 

D 

n. of duplications 


number of duplications undergone by each cell 

[0,oo] 

S 

senescence factor 


describe the senescence of each cell 

[1,0] 

Do 

n. of duplications 


n. of duplications at which senescence starts 

30 

Ds 

n. of duplications 


n. of duplications at which senescence is complete 

60 

Pn 

cell state probability 


probability that cell is in state n at a given time 

[0,1] 

Pdupl 

duplication probability 


controls probability of duplication for each cell 

[0,1] 

Parr 

arrest probability 


controls probability of arresting a cell for DSB accumulation 
or senescence 

[0,1] 

Pdeath 

death probability 


controls probability of apoptosis when DSBs in a cell approach 

N u 

[0,1] 

a 

retarding parameter 


^^cnt 

slows cell killing by DSB accumulation 

> 1 

Prestart 

restart probability 


controls probability of restarting cell cycle from an arrested 
state 

[0,1] 

P 

accelerating parameter 


accelerates cell repair capability to promote restarting 

[0.5-2.] 

^u{i) 

chemical concentration 


instantaneous concentration of species fd in cell u(i) 

[0,oo] 

m(/) 

source concentration 


constant source of species /d in cell u(i) 

[0,oo] 


diffusion time 

minutes 

inverse diffusion coefficient for species jd 

[0,00] 


B-factor concentration 


instantaneous concentration of ’bystander’ pro-mitogenic fac¬ 
tor 

constant source of ’bystander’ pro-mitogenic factor 

[0,00] 


B-factor source 


[0,1] 


B-factor diffusion coefficient 


reciprocal of 6^, the diffusion time of the ’bystander’ pro- 
mitogenic factor across the cell membrane 

[0,1] 

7 

coupling parameter 


couples concentration of B-factor to cell duplication time 

[0.1-2.] 


Table 1: List of the principal model variables. Unless differently indicated, all the variables are adimensional. 


17 




